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Abstract 

A polynomial-time method to compute the truncated theta function, 
its derivatives, and integrals is presented. The method is elementary, rig- 
orous, explicit, and suited for computer implementation. The basic idea 
is to iteratively apply Poisson summation while suitably normalizing the 
arguments of the truncated theta function in each iteration. The method 
relies on the periodicity of the complex exponential, which enables the nor- 
malization of the arguments, and the self-similarity of the Gaussian, which 
allows for repeated application of Poisson summation; in other words, the 
method relies on modular properties of the theta function. Applications 
to the numerical computation of the Riemann zeta function and to find- 
ing the number of solutions of Waring type Diophantine equations are 
presented. 

1 Introduction 

Sums of the form 

Y, <?(fc)cx P (/(fc)), f(x)zC[x] (1) 

arise in areas such as number theory, differential equations, lattice-point prob- 
lems, optics, and mathematical physics, among others. For example, one en- 
counters such sums in the context of Diophantine equations and fractional parts 
of polynomials ( |Koj 1. solutions of heat and wave equations ([M]). counting of 
integer points lying close to a curve f |Huj ) . numerical integration and quadrature 
formulas ( |Ko| ) . and motion of harmonic oscillators ( |Kaj ) . Due to the impor- 
tance such sums, there exists an abundance of methods to bound them. For 
instance, Vinogradov's [V] methods supply bounds to these sums, which, along 
with some involved sieving techniques, are used in attacking Goldbach- Waring 
type problems (see [LWY] for example). 
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Despite the substantial interest in such sums, comparatively little is known 
about how to compute them for general values of their arguments. Yet, in some 
situations, it is useful to be able to compute such sums efficiently. We later 
describe two such settings, both of which originate in number theory. 

The simplest examples of (JTJ) occur when f(x) is of degree one and g(x) 
is a polynomial. There, we obtain the geometric series and its derivatives, 
where "closed-form" formulae for their values are available. The first non-trivial 
example occurs when f{x) is a quadratic and g(k) is a polynomial. Let us 
specialize to this case. Then the sums to be evaluated can be reduced to the 
form 

1 - 

F(K,j; a, b) := k ° exp(27riafc + 2nibk 2 ) (2) 

fe=0 

In this paper, using ideas that are rooted in analysis, we prove that F(K , j; a, b) 
can be numerically computed in logarithmic time in K for any a, b, and j > 0. 
For brevity, we sometimes refer to this result as the "theta algorithm." Let us 
demonstrate the utility of our algorithm. 

For many reasons, which include numerically verifying the Ricmann Hy- 
pothesis, moment questions, and, more recently, connections to Random Matrix 
Theory models, the values of £(l/2-Hi) on finite intervals are of great interest to 
number theorists. There exist several methods to compute £(l/2 + ii) with poly- 
nomial accuracy] that is, methods to obtain the numerical values of <^(1 /2 + it) 
for t > to, to some fixed positive number, with absolute error bounded by t~ c for 
fixed c > 0. An elementary method is usually derived from the Euler-Maclaurin 
summation formula; see [Ej. It has complexity 0(t). Another method relies on 
the straightforward application of the Riemann-Siegel formula. One frequently 
used version of the formula is 

C(l/2 + it) =n^e ie wj^n- 1 / 2 exp(itlogn) \ + $(t) + 0(r 5 / 4 ) (3) 

where n\ := \_y/t/(2iv)\ , and 9(t), are certain real-valued functions that 
can be evaluated in time O(logt); see [E]. 

More recently, Odlyzko and Schonhage [OSj derived a practical algorithm to 
simultaneously compute 0{T X / 2 ) values of (,{1/2 + iT + it), where t € [0, T 1 / 2 ], 
in T 1 / 2+0 ( 1 ^ time. The Odlyzko-Schonhage algorithm did not reduce the cost 
of a single evaluation of zeta. A single evaluation still consumed Oft 1 ' 2 ) time 
using the Riemann-Siegel method. Then, Schonhage [S] lowered the cost of a 
single evaluation of zeta to t 3 / 8 +°( 1 ). Later, Heath-Brown |HB] presented ideas 
that further improved the complexity of a single evaluation to t 1 / 3 +°( 1 ). 

Our theta algorithm directly leads to another and potentially practical method 
to compute ^(1/2 + it) at a single point with polynomial accuracy in t 1 / 3 +°( 1 ) 
time. The derivation is explained in a general context in [Hj and [S] (sim- 
ilar manipulations can also be found in [T], page 99). The basic idea is to 
apply appropriate subdivisions and Taylor expansions to the main sum in the 
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Riemann-Siegel formula in order to reduce its computation to that of a sum of 
0(t 1/3 ) terms of the form F(K,j; a, b). 

As another simple and direct application of the theta algorithm, we show 
how to find the number of solutions of a Waring type Diophantine equation. 
Suppose we want to find the number of integer solutions to the system 

s s+t 

^(a.fc, + fckf) - («* fc * + = ( mod M ) ( 4 ) 

i— 1 i— s+\ 

< k\, . . . , k s+t < K, ai, , a s +t , s +t are some fixed integers 

A simple calculation reveals that the number of solutions is given by 

M-l / s \ / s+t \ 

m E II w °; °< W ^ W II °; Q ' z /^. a i/M)* (5) 

/=0 \i=l / \i=s+l / 

Since we can compute F in polynomial time, the cost of computing ^ is 
(s + t)MK ( x \ This is already significantly better than a brute-force method. 
One can use the Fast Fourier Transform to compute with sufficient accuracy 
in (s + t)MK°^ + K 3 +°( 1 ) time. But this is less efficient and it requires tem- 
porarily storing large amounts of data. In the special case M = K, one can use 
Gauss sums to solve the problem in complexity 0((s + t)M). 

To gain a better appreciation of the "generic" behavior of F(K,j; a, b) it is 
useful to consider some numerical evidence. Figure 1 depicts the real part of 
F(100,0;a,6) where either a or b runs over the interval [.110, .113]. 

9?{F(100, 0; a, .13)} - 1 SR{F(100, 0; .37, b)} - 1 




Figure 1: Two plots of SR{F(100, 0; a, 6)} - 1 for selected values. 



Since the argument of e 27riak + 27rlbk [ s more sensitive to perturbations in b than 
a, it is not surprising that F exhibits much less uniformity as b varies. Also, 
the size of 3?{F(100, 0; .37, b)} - 1 is on the order of V^OO w 14.1. This is 
consistent with the behavior of a sequence of independent random variables of 
mean zero and variance 1/2. Note however that around rationals with relatively 
low denominators the behavior is markedly different; there, sudden spikes or 
troughs may occur (e.g. F(K, 0; 1/2, 1/2) = K + 1). 
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In searching for methods to compute F(K,j;a,b), one should capitalize on 
the rich structure of the theta function. The theta function, together with 
variants, occurs frequently in number theory; it is directly related to the zeta 
function by a Mellin Transform, it has a functional equation, and other mod- 
ular properties. Thus, one anticipates that a promising method to compute a 
truncated theta function would be one that directly integrates such features in 
its design. Indeed, this viewpoint leads to a natural solution of the problem. 
Before we embark on a detailed description of the solution, we state it as 

Theorem 1.1. There exist absolute constants n,\ and k 2 such that for any 
e € (0, e -1 ), any positive integer K , any non-negative integer j, any a, b G R/Z, 
and with an underlying computational model that performs arithmetics using 
0(u 2 ) bits, where v = (j ; + 1) log(K/e), the function F(K,j;a,b) can be com- 
puted with error bounded by 0{u Kl t) using O (v K2 ) arithmetic operations. The 
big-O constants are absolute. 

See the beginning of Section 3 for the definition of an arithmetic operation. 
Also, note that a bit complexity bound follows routinely from the arithmetic 
operations bound because all the numbers that occur in our algorithm have 
0(v 2 ) bits. We also show how to compute the related sums 

K 1 

G(K, j; a, b) := 2J tj exp(27riafc + 2mbk 2 ) 

k=l 

with similar error and complexity. This is done mainly for use in the separate 
paper [H]. 

We do not try to obtain numerical values for the constants «i and K2 in 
Theorem 11.11 With some optimization, they could probably be taken to be 
around 3. Also, in a practical version of the algorithm the arithmetic can be 
performed using substantially fewer than 0(v 2 ) bits and we would likely be able 
to replace v(K,j, e) with j + log(K/e). 

Let us motivate our method to compute F(K,j;a,b) when j = 0. To this 
end, recall the following application of Poisson summation due to van der Cor- 
put; see [T], page 74, for a slightly different but equivalent version. We refer to 
this application as the van der Corput iteration, although it is not conventionally 
labelled as such. 

Theorem 1.2 (van der Corput iteration). Let f{x) be a continuous real function 
with monotonically increasing derivative in (s,t). Let f'(s) — a, f'(t) = [3. 
Then 

cxp(2vri/(fc)) = X] / ex P( 2 ni(f(x) ~ vx)) dx + Ri (6) 

s<k<t a-ri<v<f3+ri 

where < rj < 1 and Ri = 0(log(/3 - a + 2)). 

The van der Corput iteration converts a sum of t — s terms to a sum of (3 — a 
terms. In order to turn this transformation into a potentially useful computa- 
tional device, we almost surely need (3 — a < c(t — s) for some absolute constant 
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< c < 1. Moreover, we must be able to compute i?i and each of the terms in 
the sum over v, which are precisely the terms in the Poisson summation formula 
that contain critical points, using relatively few operations. Still, if c > 0, the 
length of the sum over v may be of the same order of magnitude as the length 
of the original sum. Thus, the complexity of the problem appears unchanged 
unless we can handle the sum over v efficiently. However, if we require the func- 
tion e 2m f( x } to possess some favorable Fourier transform properties that allow 
us to turn the w-terms into ones suited for yet another application of ([5]), then, 
under such hypotheses, one may hope that repeated applications of the van der 
Corput iteration are possible. If they are, one can compute the original sum 
over k in logarithmic time in K . 

Such restrictions on fix) are quite stringent, which will severely limit poten- 
tial candidates for the proposed strategy. Fortunately, the choice f(x) — ax + bx 2 
is particularly amenable to repeated applications of the van der Corput itera- 
tion. Indeed, let s = and t = K. Assume that \a] < [a + 2bK\ , which is 
frequently the case. Then, the relation © becomes 

K [_a+2bK\ K 

exp(27rmfc+27rzbfc 2 ) = / exp(27riax-f 2nibx 2 — 2nivx) dx+Ri (7) 

k=Q v=fa\ 

Let us write J7]) as Si = S r + R\ . We refer to sums of the form Si as "quadratic" 
sums. First, recall the following "self-similarity" property of the Gaussian; 

/oo 
e at "* 2 dt = V^e a ' 2/ \ aeC 
-oo 

The method that we develop works as follows. It is easily shown that we can 
always reduce to a 6 [0,1) and b € [0,1/4] in ([7]). This is important, oth- 
erwise consecutive applications of Poisson summation essentially cancel each 
other out. But with b £ [0, 1/4], S r has length < K/2. Furthermore, using the 
self-similarity of the Gaussian as well as shifting of integrals to their stationary 
phase, the sum S r is transformed to a quadratic sum of the same length plus 
a remainder i?2- Finally, Ri and i?2 are computed quickly. By iterating this 
procedure at most log 2 K times, we arrive at a sum of length O(l) that can 
be evaluated directly. In other words, most of the work of the algorithm con- 
sists of computing the aggregate of the "errors" R\ and R 2 . Note that in this 
framework, varying a corresponds to sliding the sum over v whereas varying b 
corresponds to compressing, or stretching, the sum. The latter feature largely 
accounts for applicability of the van der Corput iteration in the context of this 
paper. 

For general j > 0, the method consists of at most log 2 K iterations. Each 
iteration acts on F(K,j; a, b) in the following way 

3 

F(K,j; a,b)=^2 w ltj , a ,bF (K* bKl I; a* ab , b* ab ) + R K j, a ,b + E K ,j, a ,b (8) 

1=0 
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Let us suppress dependencies to avoid notational clutter. Then, it is shown that 
there exist absolute constants h\ and K2 such that R and the coefficients wij 
can be computed using 0(v K2 ) operations, the function Eicj, a ,b = 0(v Kl e), and 
K~i < K/2. A key point is that the tuple (K*,a*,b*) does not depend on j. 
Therefore, the number of sums F we need to compute in each iteration is always 
< To elaborate, our method acts on a sum of the form Xw=o ZijF(K, I; a, b) 
in the following way 

3 j j j 

ziF{K : I; a,b) = Y^ ™l,j,a,b F { K t,b,K, k a* a<b , &* )6 ) +^2 R K,l,a,b + ^2 E K,l,a,b 

1=0 1=0 1=0 1=0 

where 

3 

U>l,j,a,b ■= ^ Z a Wl )St afi (9) 
s=l 

and Wij y a,b are defined as in (|8|). We also show that the coefficients wij do 
not grow too rapidly with each iteration in the algorithm. In fact, in Section 
3 we show that the maximum modulus of w^j over all iterations of the algo- 
rithm is 0((j + 1)4? K 2 ). This bound is rather generous but it is sharp enough 
for purposes of our error analysis. As for the sums G{K,j;a,b), they will be 
reconstructed as short linear combinations of sums of the form F. 

The paper is organized as follows. For easy reference, Section 2 lists var- 
ious contours, integrals, and other quantities defined in the paper. Section 3 
describes the typical van der Corput iteration. Each such iteration amounts 
to evaluating certain "short" exponential integrals, which in turn are converted 
to sums. Section 4 gives pseudo-code for the algorithm. Section 5 shows how 
to compute the related sums G(K,j;a,b). Finally, Section 6 contains proofs 
of various lemmas employed in previous sections. We also decided to include 
lemmas 16.71 and 16.81 in Section 6. Although these lemmas are primarily intended 
for use in the separate paper [Hj, it seems that their natural context is here as 
they are merely specializations of the theta algorithm. 



2 Notation 

For easy reference, we list expressions defined in this paper. Define the contours 

Co := {t |0 < t < K}, Ci :={K + it\0<t<K} 

C 2 := {ie 7 ™/ 4 | < t < V2K}, C 3 = {-it\0 < t < oo} 

d = {K - it\0 < t < oo}, C 5 := {te" /4 | - oo < t < 0} 

C 6 := {ie" /4 |V2^ < t < oo}, C T := {ie" 7 " 74 | < t < VlK} 

Also, let C 8 := C 2 U C 5 U C 6 and C 9 := {t \ < t < oo}. Figure 2 depicts these 
contours. Note that dots and dashed lines determine boundaries of contours; 
the latter represents lines that continue indefinitely; we make use of the usual 
set difference operator "\" when labelling contours. 
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Define the functions 

I c (K,j;a,b) := — [ t j cxp(2niat + 2nibt 2 ) dt 
K 3 Jc 

J(K,j;m,a,b) := ^ J t j cxp (-2nat - 2mbt 2 ) 1 ^^pj^ dt 

It will be convenient to let Ic{K,j;a,b) := Ic{K,j;ia,—b). If j = 0, then 
j is omitted from the list of arguments; for example, we let Ic{K;a,b) := 
Ic{K, 0; a, b), J(K;pi, a, b) = J(K, 0;pi, a, 6), and so on. 

Let [a;J denote the largest integer less than or equal to x , \x~\ denote smallest 
integer greater than or equal to x, {x} denote x — \x\ , and logx denote log e x. 
With this notation, define 

p := \a] , q := [a + 2bK\ 
pi := q — p, uj := {a + 2bK} 
wi:=[a]-a, u(K, I, e) := (I + 1) logK/e 

We let exp(x) and e x stand for the usual exponential function (and are used 
interchangeably). Finally, we define 0° := 1 whenever it occurs. 

3 The Typical iteration for F(K, j; a,b) 

An arithmetic operation means an addition, a multiplication, an evaluation of 
the logarithm of a positive number, or an evaluation of complex exponential. 
We measure computational complexity (or time) by the number of arithmetic 
operations required. This in turn can be routinely bounded in terms of bit oper- 
ations because all the numbers that occur in our algorithm have 0(v(K,j, e) 2 ) 
bits. Frequently, we abbreviate "arithmetic operations" to simply "operations." 

For any j > and e € (0, e _1 ), we say K > is large enough if it is satisfies 
the lower bound K > A where A(K,j,e) := 1000f(AT, j, e) 6 . In particular, if K 
is large enough, then e~ K = o((e/Ar) 1000 (J +1 )). We also let i>(K, e) := v{K, 0, e) 
and A(K,e) := A(K,0,e). 

We call a real pair (a, b) normalized if (a, b) e [0,1) x [0,1/4]. The nor- 
malization is important because sums are converted to integrals via Poisson 
summation. Therefore, different choices of a or b produce different integrals. 
We remark it is mainly the normalization of quadratic argument b that truly 
matters. Normalizing a so that it is in the interval [0, 1) is not critical to what 
follows; for example, it suffices to take a £ [— m, m] for some < m = 0(1). 

Let j be a non-negative integer, e any number in the interval (0,e _1 ), K 
any large enough integer, and (a, b) any normalized pair. Then, with p and q 
defined as in Section 2, either q > p or q < p. The first possibility is the main 
case, and it is where the method typically spends most of its time. The second 
possibility is a boundary point. We show how to handle each case separately. 
Arithmetic is performed using 0(y(K, j, e) 2 )-bits. 
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3.1 Case: q > p 

By Poisson summation we have 

F(K,j;a,b) = i (5j + exp{2niaK + 2iribK 2 )) + PV ^ I Co (K,j;a~m,b) 

m=—oo 

where Sj is Kronecker's delta, and PV stands for Principal Value. Define 

i 

Si(K,j;a,b) := ^ I Co (K,j;a - m,b) 

m—p 

S 2 (K,j;a,b)-=PV £ I Co (K,j;a - m,b) 

By Lemma 16.11 the condition q > p implies 2bK > 1. This simple observation 
is used repeatedly throughout Section 3.1, often without any comment. 

3.1.1 The sum Si{K,j;a,b) 
By Cauchy's Theorem 

Ic (K, j;a-m,b) = Ic 2 (K, j; a - m, b) - I Cl (K, j; a - m, b) (10) 

Consider the integral over C\ first. Ignoring the term m = q for now and 
summing over m £ {p, . . . , q — 1} we obtain 

^2 I cA K -,3\a-m,b) = ciJ^iM J J(iT,Z;pi,w,6) 

m=p 1=0 ^ ' 

where c\ = i exp(2iTiaK + 2iribK 2 ). When m = q, Cauchy's Theorem yields 
I Cl (K,j;a-q,b) = Cl Yy{^jIc a {KJ^M 



Cl 



i l [\ )lc 7 {K, /; w, 6) - Cl ]T *' ( \ ) *5 &) 

Z=0 ^ ' i=0 



According to lemmas 16.21 and 16 . 31 there exist absolute constants £3 and £4 such 
that for any e £ (0,e _1 ), any integer I > 0, any integer M — 0(e 2K ), and any 
integer K > A(K, l,e), the functions J(K,l; M,u,b) and I 7 (K,l;uj,b), where 
7 £ {C7, Ci}, can be computed in 0(v{K, j, e) K3 ) operations with error bounded 
by 0{v{K,j,t) k ^-iK- 2 e). 

Having disposed of the integral over C\ in (|10p . we now consider the one 
over C2. Write 

Ic 2 [K, j',a- m, b) = I Cs {K, j;a-m,b)- I Cs (K, j; a - m, b) - I Ce (K, j;a-m,b) 
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Consider C5 first. Cauchy's theorem and a straightforward estimate yield 
8 

^2 Ic 5 (K,j;a-m,b) = c 2 J (K,j;p 1 ,uj 1 ,b) + 0{e~ K ) 

m=p-\-l 

where C2 — (— l) J e^ +1 ^ 7r ^ 2 . When m = p we have 

7 C5 (X, j;a-p,b) = cjc 7 (K, j> x , 6) + O(e ^) 

As before, the integrals J(.) and /(.) in the last two equalities are handled by 
lemmas [6.21 and 16731 As for Cq, it is not hard to see that Ic 6 (K,j;a — m,b) is 
0(eT K /K) for each m = p, . . . , q — 1; hence, these terms are negligible. When 
r7i = q, we obtain 

Ic 6 (K, j;a-q,b) = c 3 e- 2 ™ K ^ Q 2^/c 9 Z; w - 7^ + 26^ + 726^, -2*6) 

where c 3 = e 0'+i )™/4+2"aK _ xhe i nte grals J Cs) (K,l;ui-iuj + 2bK + i2bK, -2*6) 
are handled by Lemma 16.31 Finally, by Lemma l6.4l 

s j i 

X! I c a {K,j\a~ 777,6) = ^ui S!i F((7,s; «!,&!) - g x _ p io aj - (11) 

m—p s—0 s—0 

where ai = a/(26)(mod 1), 61 = — 1/(46) (mod 1), and all w s j can be computed 
using 0(j A ) operations. 

For purposes of our error analysis, it is sufficient to find the maximum mod- 
ulus of w St j during a full run of the algorithm. Suppose 2bK > A(K,j, e), then 
by Lemma 16.51 

J_ e l+jyA gl+l/ log 2 K 

\w., m \ < — >=-(! + 2 i/A) < 7=— (1 + 1/loga A - ) 

V26 V26 

Now, suppose < 2bK < A(K,j,e). This can happen only in the last iter- 
ation. There are two possibilities: either p < q < A(if, j, e) + 1 or q < p. In 
the latter case we do not reach (fTTj) at all. In the former case, Lemma 16.51 
gives YL= S K,m| < (j + l)4^+ 2 (26)- 1 / 2 . Recall ©. Then, put together, the 
maximum modulus of the coefficients w s j that can occur over a full run of the 
algorithm satisfies the bound 0((j + l)4^Vz?e log 2 K ) = 0((j + 1)4-? if 2 ). This 
concludes our computation of S\(K,j; a, 6). 

3.1.2 The sum S 2 (K, j; a, 6) 

Let us first handle the subsum X)m=<j+i ^c (K, j; a—m, 6), M > q. For m > q+1 

\I Co -iT(K,j;a-m,b)\ < (2Ty e - 2 ^'^ T [ e-^^dt -r^, (12) 

Jo 
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So, Cauchy's Theorem gives 

Ic (K, j;a~m,b) = Ic 3 (K, j; a - to, b) - I Ci (K, j; a - m, b) (13) 

We remark if j = 0, then (fT2"|) holds uniformly in a € [0, 2] and to > q. Therefore, 
(fT3"|) holds for all a £ [0, 2] and to > q. This observation is used the proof of 
Lemma 16.71 By a simple calculation 

M 

J2 Ic 3 (K,j;a- to, 6) = Ci J(K,j;M - q,2bK - u,b) + 0(e~ K ) 

m= g+1 

where C4 = (— i)' J+1 . A similar calculation gives 

m j , . 

J2 IcAKJ; a-m,b) = c 5 ^H)' ( j ) W l;M - q-1,1- u,b) + 0( e - K ) 

m=q+2 1=0 ^ ' 



with c 5 = _j e 2™»K'+2iri&ir _ Furthermore 

IcAK, j;a-q-l,b) = c 5 ^(-z)' Q J Cg (X, I; 1 - w , 6) 

The integrals Ic 9 {K, I; 1 — w, 6) are handled by Lemma HT51 

As for Ylm=-M Ic (Ki 3'i a ~ m ib), the situation is analogous. Simply use 
the conjugates of the contours C3 and C4, then repeat the same calculations 
with the appropriate modifications. The resulting integrals are 

p-2 

^2 I c^( K 'j' a ^ m ' b ) =C6J{K,j;M + p- l,l-wi,6) 

m=—M 

2 i; o - m, 6) = or 53 ( j ) W M + P. 2Wf ~ w i > 5 ) 

%i"(-^> J! a ~ P + & ) = c 6^c 9 (K, j; 1 - Wi, b) 

where cq := i- 7 ^ 1 and C7 := j e 27rm - ft: + 27rlf ' i <' , Finally, the terms corresponding to 
\m\ > M can be bounded as follows. Write 

2 [ K 

PV 53 Icq(Ki 3\ a ~ m i b) — 53 I e-xjp(2iriax + 2iribx 2 ) cos(2irmx) dx 

\m\>M m>M ^° 

Integrating by parts this is equal to 

i-K 



[ — - — 7 / (1 — 5j)ocP 1 exp(27riax + 2itibx 2 ) sm(27rmx) dx+ 

£it\™ K3 Jo 

J x 1 (a + 2bx) eyLp(2iriax + 2nibx 2 ) sin(27TTO2;) dx j 



m>M 

2z ' A ' 
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By the Second Mean Value Theorem, we deduce for M > 2K 



Py Ic (K,j;a-m,b) = olYl 

\m\>M \rn>M 



m(m — K) J \M 
Finally, take M = e 2K to obtain the bound 0(8~ j K~ 2 e~ K ) . 
3.2 Case: q < p 

The case q < p is a terminal point of the method. Observe if q < p, then 
< a + 2bK < 2. Thus, we may assume l/K 2 < b < 1/K (if b < 1/K 2 , then 
the problem is trivial). We may also assume K is a multiple of 8. Write 

e 2iriaK+27ribK 2 Y 7 (m+l)K/8-l 

F(K, j; a, b) = — + — ^ £ k J ex V (2niak + 2mbk 2 ) 

m=0 k=mK/8 

But 

(m+l)K/8 j . . 

— Y fc J exp(27rwfc + 27ri6fc 2 ) = c 13 8~ j ^2m j - l ( j \ F(Ki,l;a,b) 

k=rnK/S l=0 ^ ' 

where | C13 1 = 1, < m < 8, K\ := K/8, and a := a + mbK/A. We can normalize 
so that -1/2 < a < 1/2. Since < 2bK < 2, then < 2bK 1 < 1/4. So, 
< \a\ +2bKi < 3/4. Put together, we may now assume that \a\ + \2bK\ < 3/4, 
where a, b £ R. Define the function 

x a 

f a (x) := exp(2niax + 2wibx ) 



where a is a non- negative integer. By Lemma 16.61 

max \fi N \x)\ < ( + 2n(\a\ + \2bK\) ^ 

0<x<K \ A 

So, given a sum F(K,j;a,b) and e G (0,e _1 ) such that K > A(K 7 j,e), a, b e R 
and \a\ + \2bK\ < 3/4, we can apply Euler-Maclaurin Summation to F(K, j; a, b) 
with N correction terms. Take N = log(8 J A' 3 /e)/(2 log(8/7)). Then, resulting 
absolute error is bounded by 

K 

\f\ 2N+1 \x)\dx < 2K(7/8)~ 2N = 0(8- j K~ 2 e) 



{2n) 2N 

It remains to evaluate the integral Iq (K, j; a, b). But this is handled by Lemma [6~3l 
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4 The Algorithm for F(K, j; a, 6) 

Lemma 4.1. For any integer K > 0, any integer j > 0, and any a, o G C, i/ie 
junction F(K, j; a,b) satisfies the identities 

F(K, j; a, 6) = F(K,j; a + 1, 6) - j; a, 6 + 1) = F(K,j; a ± 1/2, 6 ± 1/2) 

/Yoo/. This follows from e 27r4 ^ +1 ) = e 2T " z and (fc 2 ± fc)/2 € Z for fc € Z. □ 

Lemma 14.11 implies that for any integers K > and j > 0, there exists a pair 
(ao,&o) G [0,1) x [0,1/4] such that F(K,j;a,b) is equal to F(K,j;a ,bo) or 
its conjugate (this is independent of K and j). We are now ready to present 
pseudo-code to compute J2l=o z i,jF{Ki h a , &)i z ij — 0(1). It suffices to do the 
arithmetic using 0{v(K, j, e) 2 ) bits. 

INPUT: a, 6 e M/Z, an integer K > 0, a positive number e < e _1 , an integer 
j > 0, and an array of complex numbers zi t j, I = 0, . . . ,j such that 
zij = 0(l). 

OUTPUT: a complex number sum that equals X)z=o z i,jF(K, I] a, b) modulo 
error 0(v(K,j,e) kl+1 e), where k\ is an absolute constant. 

INITIALIZE: sum = 0, flag = 0, counter = 

1. Normalize (a, b) — > (ao, £»o)- Set flag = 1, zjj — > z; j if conjugation is 
needed to normalize (a, b) (see Lemma T4. II) . 

2. Set p = \a ] , g = [a + 2b K\ . 

3. If K is large enough, then go to 4. Otherwise, compute the sum directly, 
let R[counter] be the result. If flag = 1, set R[counter] — > R[counter\. 
Go to 9. 

4. If g < p, then apply the typical iteration in this case, let R[counter] 
denote the result. If flag = 1, set R[counter] — » R[counter\. Go to 9. 

5. Apply the typical iteration in the case g > p 

3 

F{K,j; a , bo) -> J2 v>ljF{Ki,l\ °>i, M + Rj + 0{v{K, j, ef'e) 

1=0 

i-i 

F{K,j - V,aoM) ^^mj-xFiKx^aiM) + Rj-i + 0{v{K,j, 

1=0 

F(K, 0; a , b ) -» iflo.o^i, Z| ox, &i) + i?o + 0(v(K,j, e)~ Kl e) 
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6. Set R[counter] = J2i=o z i,j R ^ z hi ~* YH=i z s,] w i,s, 

K — ► K\ = [ao + 2boK\ , counter — > counter + 1, a — ► Oi, b — ► 61. 



7. If /^aa = 1, then set zjj — > R[counter] — > ii [counter], a — > —a, 

— > —6, //ay — * 0. 

8. Goto 1. 

. sum — > 2^/=o ^IAr 

5 The Sums G(K,j;a,b) 

We show how evaluate the sums G(K,j; a, b). Let TV be a positive multiple of 
16. Define 

2JV-1 

V(JV, j; a, 6) := ^ -j exp(2?riafc + 2nibk 2 ) 

k=N 

Since G(K 7 j;a,b) consists of O(logi^) sums of the form V(N,j;a,b) plus a 
remainder sum of length 0(logK), it is enough to compute V. We can write 

15 JV+(m+l)A7l6-l 

V(JNT,j';a,6) = ^ ^ — exp{2niak + 2nibk 2 ) 

m=0 fe=W+mAT/16 



— exp(27riafc + 27ri6fc 2 ) (14) 



Define N m := N + mN/W. Then 

N m+1 -1 

h' 

k=N m 

00 /• , 7_1\ N / 16 ~ 1 U 

= ^ L E(- 1 )'rt-T ) £ ^ r ex P (2 7 r^ + 2rf) 

where |ci. m | = 1. Since ( 3 ^_l )k l /N^ < 8~', we can truncate the sum over I 
after v(K,j,e). Finally, the inner sums in l|14p are handled by Theorcm ll.il 

6 Auxiliary Results 

Lemma 6.1. Let a, b £ R and K > 0. If [a + 2bK\ > \a\ then 2bK > 1. 
Proof, [a + 2bK\ > \a] => 2bK + a > \a\ + 1. So, 2bK > 1. □ 

Lemma 6.2. Let J(.), f(-), and A(.) 6e defined as in Section 2. There exist 
absolute constants K3 and K4 such for any e € (0, e _1 ), any integer I > 0, any 
positive integer K satisfying K > h.(K, l,e), any b € [0, 1] satisfying 2bK > 1, 
any < cj < 1, any w > 0, and any positive integer M satisfying M — 0(e 2K ), 
the integral J(K, I] M, w + to, b) can be evaluated with error bounded by 
0(i'(K,l,e) K3 8~ l K~ 2 e) using 0(v{K,l,e) Ki ) operations. Big-0 constants are 
absolute and arithmetic is done using 0(v(K, I, e) 2 ) bits. 
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Proof. Let us assume w = 0, the case w > is similar. Let L = \3i/(K, I, e)] . 
Then 



J(K,l;M,cj,b) = Iff t l exp (-2mjt - 2mbt 2 ) 1 GXP ( ^ Mt) dt 
K l J v ' exp(27rf) — 1 

+0(8^A^ 3 e) 

So, it suffices to evaluate the integrals 

1 f l+1 , i „ „ , 2N 1 - exp (-2?rMi) , 

-tt / r exp ( -2?rwt - 27ri6t 2 ) P— dt 

K l J n t V 7 cxp(27rf) - 1 

where n = 0, . . . , L — 1. Taylor expansions yield 

exp(— 27rijjn — litibn 2 ) s—^ ( l\ l _ s ^-y (— 2-n:ib) r 



dt + ((log M)8~ l K 



s=0 v 7 r=0 

r^+ar / o . a ., .x l-exp(-27rM(t + n)) 

/ i s+ ^ r exp (-27tcj£ - Anibnt) f-^- ^ — 

Jo exp(27r(t + n))-l 



In order to evaluate the last expression, it suffices to deal with the integrals 

1 n. / « , , s 1 - exp(-27rM(t + n)) , , . 

t a exp (-2nut - 4nibnt) ±-- £ ^ dt (15) 

o i V ; exp(27r(i + n)) - 1 V 7 

for integers a S [0, 3i]. These can be evaluated by unfolding the geometric 
series; that is, write (fT5|) as 



M 

r(m + u + 2ibn)t) dt (16) 



M 1 

exp(— 2irmn) / i" exp (— 27r( 
m=i • /o 



Apply Taylor expansions to the factor exp(— 2iruit) in (JTHJ) to reduce its 
evaluation to that of sums of the form 

M .1 

^ exp(-27rmn) / r +/3 exp (-2w(m + 2ibn)t) dt (17) 
m=i -^o 

for integers < f3 < L. Let a.\ := a + /?. For m > «i + 1, we explicitly 
evaluate (jTTJ) to obtain 



E +1 , .,,,1 ai! \— * . exp(— 27rm — 47r6in)) 

v ; (ai + l-w)! ^ y ' (2mn + Anibn ) v 

v=l v y m=ai+l v ' 

,/ lV n i exp(-27rmn) 

+ ai! 2. (2. TO + 4«6n)« 1+ i (18) 

m— ai + l 7 



The sums in (jT8|) can be evaluated by truncation or Euler-Maclaurin 
summation (in the case n = in the second sum). Now, consider the terms 
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1 < m < a\. For each integer < n < L and each integer 1 < m < ct\, there 
are two possibilities: either n + l<m<aiorl<m<n+l. In the first 
case, apply the change of variable t — > mt to the integral in (fT7|) to reduce it to 
integrals of the form 



1 



t ai exp (-2vr(l + 2ibn/m)t) dt 

where / is an integer in [0,m]. These integrals are straightforward to evaluate; 
first, make a change of variable i — s- 1 — Z, then use Taylor expansions to reduce 
the integrand to a polynomial in t of degree O(L). Finally, the case 
1 < m < n + 1 is handled analogously. □ 

Lemma 6.3. Let Ic(-)> Ic{-)> v{-)> C\, Cj, and Cg be defined as in 

Section 2. There exist absolute constants and Ke such for any e G (0, e _1 ), 
any integer j > 0, any positive integer K satisfying K > A(K,j, e), any 
b € [0, 1] satisfying 2bK > I, any real a satisfying a = 0(1), any < u> < 1, 
and any < 5 < 1, each of the integrals 

Icr(K,j;uJ,b), Ic T (K,j;ui,b) 
I Cg (K, j; 5 + uj,b), I Cg (K, j; uj - iw + 2bK + i2bK, -2ib) 

can be evaluated with error bounded by 0(v(K, j, e) K5 8~ l K~ 2 e) using 
0(y(K, J, e) Ka ) operations. Under the same assumptions, except now b need 
not satisfy 2bK > 1, the integral Ic (K,j;a,b) can be evaluated with the same 
error and complexity. Big-O constants are absolute and arithmetic is done 
using 0(v(K,l,e) 2 ) bits. 



Proof. We evaluate I-^(K,j;uj,b) first. We have 

I—{K,j; uj, b) = c 8 e- 2 ™ K V ^j- [ t l exp (2mut - AirbKt + 2nibt 2 ) dt 

i=o W K J o 

where c% :— — i exp(— 2mbK 2 ). By truncating at L — \3v(K, j, e)] , the 
evaluation of I-Q-{K,j; to, b) is reduced to that of integrals of the form 

1 f n+1 

-j / V exp (2niujt - AirbKt + 2nibt 2 ) dt (19) 

for integers < / < j and < n < L — 1. In order to compute (fT9|l . substitute 
t — > t — n, then eliminate the quadratic term 2iribt 2 using Taylor expansion. 
This results in easily-calculable integrals of the form 

/ t a exp (2Tri(uj + 2bn)t - AnbKt) dt 
Jo 

where a is an integer in [0,3L]. The evaluation of 

Ic 9 (K, l;uj-iu + 2bK + i2bK, -2ib) 



1G 



is similar to I^(K,j;u),b). So, we may move on to Ic 7 (K, j;uj,b). By 
definition 



Ic T (K,j;u,,b) = H 



f V2K 

J t 3 exp {-\f2-Kuot + ^/2mujt - 2irbt 2 j dt 



where eg = exp (— (j + The change of variable t \J~bt yields 

So, truncating the integral at \VL~\ reduces the problem to evaluating 



1 



n+l 



t 3 exp -2-K—==t + 2m— =t - 2-Kt 2 dt 



mWKJj n "V W2b V2b 

for integers < n < \\/L~\ . Finally, these integrals are handled as before; 
substitute t —* t — n, then eliminate the quadratic term using Taylor 
expansions. This procedure results in readily-calculable integrals. Next, we 
consider Ic g (K,j;5 + iv,b). Let ui' := ui + 5. Then < u>' < 2 and 



1 



T 

(T - it) 3 ' exp (-2to'(T - it) - 2mb(T - it) 2 ) dt 



(20) 



So, we can replace Cg by e m ^Cg in Ic 9 (K,j;uj',b). A simple estimate then 
yields 

I e -^ / ±c 9 {K,3\u J \b)=I C7 {K,3-u l ,b) + 0{e- K ) (21) 

We have already shown how to compute the right side of (f2"Tj) . We remark if 
j = 0, then (|2H)) . hence (|21[) . holds uniformly for ui' 6 [0, 1]. This observation is 
used in the proof of Lemma 16.71 

Finally, we consider the integral Ic (K,j; a,b). This may contain a 
critical point or it may not according to whether —a/ (2b) £ [0, K] or not. We 
showed how to deal with these possibilities in Sections 3.1.1 and 3.1.2 
respectively. □ 

Lemma 6.4. Let Ic(-) an d Cg be defined as in Section 2. For any integer 
K > 0, any integer j > 0, any integer m, any a£l, and any b > such that 
q := [a + 2bK\ is not zero, we have 



Ic 8 {K,j;a- m,b) = exp ( ^-m- ^f-m 2 J ^ 



V 2b 46 ^ q 

' s=0 H 



J 



(22) 



wher 



Ws ' 3 ~ q Z>/ 2 s\(2y/bn) i+1 Ki IV T 



i 



^4 %- s -Q( mo d2)(-l) 0+ ' s)/2 / 3 „ /4 /2tt i 
^ 1 V b I 1 dJ 

i=0 '* 2 • \ / 
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We remark that \22\l is what one would expect; it is also essentially 
independent of K . 

Proof. By a change of variable 

e 0'+l)wi/4 pea 

Ic a (K, j; a — m,b) — — — — : / t-' exp(iat — at — bt 2 /it — irat + rat) dt 



(V2tt)J+ 1 ^ J-c 
e 0+i)W4 r°o (i-J%a b 



V^F(m - a)(l - i) _ V^e- 3 "/ 4 (m - a) 
26 ~ Vb 



t 1 exp t 1 1 dt (24) 



where 



Let Hj(x) denote the Hermite polynomials, which are defined by generating 
function 

di 

H^a) := (-iy exp (a 2 /2) exp (-a 2 /2) 

By the absolute convergence of (f2"4")l , it is possible to differentiate under the 
integral sign. Thus, we obtain 

e {j+l)ni/4 f d J f°° fi^2b , , 2 



Ia,(K,j;a — m,b) = — = — = - — r / exp — —at — bt /n 

8V ; (V^iry+iRj \iV2b) daJ - 



o 7ri/4+37rij/4 

exp (-a 2 /2) Hj (a) 



{2Vbny+ 1 Ki 

The coefficient of x l in Hj(x) can be interpreted as (— 1)0'— 0/ 2 times the 
number of unordered partitions of a j-element set into Z singletons and 
(j — Z)/2 unordered pairs (see [I]). Therefore, 



IcAK, j: a — m,b) — = exp — m 77-m x 

The sum over Z can be rewritten as 

Finally, change the order of summation in ([25]) to obtain the result. □ 

Lemma 6.5. Let A(.) and v(.) be defined as in Section 2. For any 

e G (0, e _1 ), any a G [0, 1], any b G [0, 1], any integer j > 0, any positive 
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integer K satisfying K > A(K,j, e), any integer < s < j, and with w s , m 
defined as in A23\) and under the condition \a\ < [a + 2bK\ , we have the bound 

m— s v \ / g—0 v y 

Moreover, if 2bK < v{K,j,e), then YL=s \ w s t m\ < (j + l)^ j+1 {2b)- 1 / 2 . 
Proof. To obtain the first bound note that 

< I s ^ (m+£T y*. (v^F|a|)'6(« 



(2wo-^6^ (^ m (2^- £T ^ J " ^ 



V26V 2bK J ^A?bK 

v n / m=0 v 

If 2M£T < v(K,j, e), then since A" > A(K,j, e), it follows that 6 < 1/(2 j + 2) 2 . 
Therefore 

f I | < 2g s ^ (m + s)! (j + 1)4^ 

~ - (2feA)- % /2& j ^ s!m!(26A)™ - ^2& 

□ 

Lemma 6.6. For any integer a > 0, any integer m > 0, any integer K > 0, 
and any rea? numbers a and b, the function 

x a 

f a (x) := — — exp(27riax + 27ri6x 2 ) 

satisfies 

max < (27r(|a| + \2bK\) + (m + a)/A) m 

0<:r<A 

Proof. Since a, 6, and AT are fixed, we suppress dependencies on them. We can 
write 

f£ n \x) = P m (x) exp(2max + 2mbx 2 ) 

where 

m+a 

J 



P m {x) := di 



m,aX 




for some complex coefficients d;, m , a . Define |P m (a;)|i := J^ST* Mz.m.a^ 
Suppose 

max |P OT (a:)|i < (2tt(|o| + \2bK\) + (m + a)/A) m 

0<IE<JC 
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Then 



max |P m+1 (a;)|i < max \2niaP m (x)\ 1 + max \AmbxP m (x)\\ 

0<x<K 0<x<K 0<x<K 

+ max |P^(a;)|i 

0<x<K 



< {2n{\a\ + \2bK\) + (m+l + a)/K) 



Lemma 6.7. For any a £R, any b > 0, and any K > let 



m+l 



□ 



q a = [a + 2bK\ , p a = \a] , p 1>a = q a - p a 
Lu a = {a + 2bK}, uj ha =p a -a, M = e 2K (26) 



Cl.a — S 2 - Pa , C2, a = 5 qa -K 1 -\, Cj,.a — 5q a -Ki-2 (27) 

f(x) = e 2 ™ K , g(x) = f(x)e- 2 * K "*, h(x) = -±= e -^IW) 

Also, let F(K;a,b) := F(K,0;a,b), K* = \2bK\, a* = a/ (2b), and 
b* = —1/(46). Then for any K > 1000 and any tuple (a, a, b) in 
[0,1] x [0,2] x [0,1/4] satisfying 

Pa+ax < q a +ax, for all X G [-1/4, 1/4] 

(a + ax,b) e (0,2) x [0,1/4], for all x e [-1/4, 1/4] (28) 

we have 

F(K;a + ax,b) = e™ /4 h(a + ax)F (k* ;a* + ||,&*) + R(K, a + ax,b) 
+0(e- K ), for x e [-1/4,1/4], 

where R(K, a + ax, b) is a linear combination of the functions (defined in 
Section 2) 

J (K; Pl.a+ax , Ul,a+ax , b), f(ax)J(K;pi, a+ax , LUa+ax, V) 

J(K; M, 2bK - uj a+ax , b), f(ax)J{K; M, 1 - w a+ax , b) 
J(K;M,l-oj^ a+ax ,b), f(ax)J{K;M,2bK-u ha+ax ,b) (29) 

as well as the functions (defined in Section 2) 

g(a + ax)I Co {K; -iiv a+ax + 2bK, -b), g(a + ax)I Co (K; e m/i (-iu) a+ax + 2bK), 
Ic 7 {K; u>i, a+ax , b), f(ax)Ic 7 {K; uj a+ax , b) 

b), f(ax)I Cr (K;l-uj a+ax ,b) 
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and the functions 

C2, a+ax h(* + ax)e 2 ^ K ^^ 2b \ c 3 , a+ax h(a + ax)e 2 * ia ^ +1 VW 

c*, a+ax h{a + ax)e 2 ^ K ^ 2b \ c ha+ax h(a + ax)e 2 ™*'^ 

f{ax), h(a + ax) (31) 

as well as the constant function. The coefficients in the linear combination can 
all be computed using 0(1) operations, are bounded by O(l), and do not 
depend on x. Big-O constants are absolute. 

Proof. This follows from our main result and the remarks following formulas 
<P1) and dHJ). □ 

Lemma 6.8. Let v(K, e) and A{K, e) be defined as at the beginning of Section 
3. Recall the definitions i26\) . $27\) , and the definition K* :— [2bK\ in 
Lemma \6.7\ Then, for any e € (0, e _1 ), any K large enough, any interval 
(w,z) C [—1/4,1/4] and tuple 

(a, a, b) e [0, l/A(K, e)] x (0, 2) x [0, 1/4] 

satisfying assumptions H28\) as well as the assumption 

Pa+ax and q a +ax are constant over x £ (u>, z), 

and for x € (w, z), there exist constants Xg , Ag € [0, 1], 9 £ {—1, 1}, 
independent of x and satisfying < A^ + 6ax < 1 for x <E (w, z), such that 
each of the functions in 129\) . 130\) . and 131\) is equal to a linear combination 
of functions of the form 

m m 2-iriaKx 2-Ki^M-x — 2iri^r-x 2 

(A ± + ax ) m e 2mL ( x t+ 0ax )-^ R ( x t+ eax ) i 

e 27T2N(\f+e a x)-27T(l-i)m / t ; e -27r(l-i) » -j= t-2nmt ^ 

JO 

plus an error bounded by 0{A(K,e)K~ 2 e). Here, < I, m — 0{v{K,e)) are 
integers and 

Ne{0,K}, M e {-1,0, K*,K* + 1}, 
K<L<K + u, K < R < K + v 
< u = 0(y{K, e)) < v = 0{v{K, e)) 

The length of the linear combination is 0{A{K,e)). The coefficients in the 
linear combinations can all be computed in 0{A{K, e)) operations, are bounded 
by 0(K), and are independent of x. Big-O constants are absolute and 
arithmetic is done using 0(v(K,t) 2 ) bits. 

Proof. This is readily deducible from the proofs of lemma 16. 21 lemma 16.31 and 
the assumption that p a +ax and q a + ax are constant over (w, z). □ 
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